package com.numericalmethod.algoquant.model.ralph2009.draft;

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.text.SimpleDateFormat;
import java.util.Date;

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;

import com.numericalmethod.algoquant.model.ralph2009.Brandt2005Regression;

/**
 * Trial implementation of 
 * 
 * "A Simulation Approach to Dynamic Portfolio
 * Choice with an Applicaion to Learning About Return Predictability"
 * (Brandt et al, 2005)
 * 
 * with one risky asset and no state variable.
 * 
 */
public class OneAssetSimulationTrial2 {
	
	
	// helper functions
	private static double average(double[] a) {
		double result = 0.0;
		for (int i=0;i<a.length;i++) 
			result+=a[i];
		result /= ((double) a.length);
		return result;
	}
	
	private static double stdev(double[] a) {
		double avg = average(a);
		double sum = 0.0;
		for (int i=0;i<a.length;i++) 
			sum+=Math.pow(a[i]-avg, 2.0);
		sum /= ((double) (a.length-1));
		return Math.sqrt(sum);
	}
	
	// CRRA Utility
	private static double utilityFunction(double w, double gamma) {
		
		return Math.pow(w, 1.0-gamma) / (1.0-gamma);
	}
	
	// (TESTING DISPLAY)
	private static void printComparison(double[] nextExcessReturn, double estimatedOptimalWeight, 
			double mu, double sigma, double riskFree, double dt, double gamma) {

		System.out.println("==== COMPARISON ====");

		double averageExcessReturn = average(nextExcessReturn);
		double returnVolatility = stdev(nextExcessReturn);
		System.out.println("Next Step Excess Return = "+averageExcessReturn);
		System.out.println("Next Step Vol = "+returnVolatility);
		System.out.println();
		
		// close form for Mean/Variance portfolio optimization on CRRA Utility 
		// ???CORRECT????
		double closeFormWeight = 1.0/gamma * ((mu-riskFree)/(sigma*sigma));
		
		// (TESTING DISPLAY)
		// the portfolio return mean and vol using our estimated weights
		double ptfMean = averageExcessReturn*estimatedOptimalWeight+(1.0-estimatedOptimalWeight)*riskFree*dt;
		double ptfVol = returnVolatility * estimatedOptimalWeight;
		System.out.println("Using Simulation Estimated Weight: "+estimatedOptimalWeight);
		System.out.println("portfolio mean = "+ptfMean);
		System.out.println("portfolio vol = "+ptfVol);
		System.out.println("portfolio mean/variance = "+(ptfMean/(ptfVol*ptfVol)));
		System.out.println();

		// the portfolio return mean and vol using the close form weights.
		ptfMean = averageExcessReturn * closeFormWeight+(1.0-closeFormWeight)*riskFree*dt;
		ptfVol = returnVolatility * closeFormWeight;
		System.out.println("Using Mean/Variance Analytic Weight: "+closeFormWeight);
		System.out.println("portfolio mean = "+ptfMean);
		System.out.println("portfolio vol = "+ptfVol);
		System.out.println("portfolio mean/variance = "+(ptfMean/(ptfVol*ptfVol)));
		System.out.println("====================");
		System.out.println();
		
	}
	
	private double[] optimalWeight;
	
	public OneAssetSimulationTrial2() {
	}
	
	public double[] getOptimalWeight() {
		return optimalWeight;
	}
		
	
	/**
	 * Solve for optiomal weights
	 * 
	 * @param rand random generator to be used
	 * @param mu annualized return rate
	 * @param sigma annualized volatility
	 * @param riskFree risk free rate
	 * @param pathCounts simulation counts
	 * @param stepsCount the algorithm requires at least 2 steps
	 * @param gamma CRRA utility parameter
	 * @param dt time step size
	 */
	public void solve(final double gamma, final PriceDynamics priceDynamics) {			
		
		optimalWeight = new double [priceDynamics.timePeriod-1];
		double[] optimalWeight2 = new double [priceDynamics.timePeriod-1];
		
		double[] aTPlus1Vec = new double[priceDynamics.pathCount];
		double[] bTPlus1Vec = new double[priceDynamics.pathCount];
		double[] aMyopicTPlus1Vec = new double[priceDynamics.pathCount];
		double[] bMyopicTPlus1Vec = new double[priceDynamics.pathCount];		
		double[] vtVec = new double[priceDynamics.pathCount];

		double[] regATPlus1Vec = null;
		double[] regBTPlus1Vec = null;
		double[] regAMyopicTPlus1Vec = null;
		double[] regBMyopicTPlus1Vec = null;		
		
		//initialization for t=T-1
		RealVector factor1Vec = new ArrayRealVector(priceDynamics.pathCount); 
		factor1Vec.set(1);//set to 1 for t=T-1
		double[] factorVec = factor1Vec.toArray();

		double[] averageExcessReturn = new double [priceDynamics.timePeriod-1];
		double[] excessReturnVol = new double [priceDynamics.timePeriod-1];
		double[] averageExcessReturn2 = new double [priceDynamics.timePeriod-1];
		double[] excessReturnVol2 = new double [priceDynamics.timePeriod-1];

		double[][] excessReturnMatrix = priceDynamics.excessReturnMatrix.getData();
		double[] phiTPlus1  = null;
		for(int t=priceDynamics.timePeriod-2; t>=0; t--){
			double[] excessReturnVector = priceDynamics.excessReturnMatrix.getColumnVector(t).toArray();
			double[] divYieldVector = priceDynamics.divYieldMatrix.getColumnVector(t+1).toArray();
			
			//determine function phi(t+1) and stored in factorVec using eqt.5 p.7 in brandt(2005)
			for(int i=0;i<priceDynamics.pathCount;i++){
				if(t < priceDynamics.timePeriod-2){
					double factor = (optimalWeight[t+1] * excessReturnMatrix[i][t+1]) + (1+(priceDynamics.riskFree*priceDynamics.dt));
					factor = Math.pow(factor,1-gamma);
					factorVec[i] *= factor;
				}
			}
			
			
			phiTPlus1 =  Brandt2005Regression.estimateYValues(excessReturnVector, divYieldVector, factorVec);
			factorVec = phiTPlus1;
			
			for(int i=0;i<priceDynamics.pathCount;i++){
				//Eq.12, redefine At+1 = E[phi(Zt+1)*Ret+1] and Bt+1 =E[phi(Zt+1)(Ret+1*Ret+1)]
				//This is the case when innovation of state variable Zt+1 coupled to the optimal weight determination
				aTPlus1Vec[i] = factorVec[i]*excessReturnVector[i]; 
				bTPlus1Vec[i] = factorVec[i]*(excessReturnVector[i]*excessReturnVector[i]);
				
				//if Zt+1 is independent of Ret+1 => reduce to myopic strategy
				aMyopicTPlus1Vec[i] = excessReturnVector[i]; 
				bMyopicTPlus1Vec[i] = (excessReturnVector[i]*excessReturnVector[i]);
			}
			
			
			//refit At+1 and Bt+1 for optimal weight determination.
			System.out.println("Fitting At+1 at t="+t + "...");
			regATPlus1Vec = Brandt2005Regression.estimateYValues(excessReturnVector, divYieldVector, aTPlus1Vec);
			regAMyopicTPlus1Vec = Brandt2005Regression.estimateYValues(excessReturnVector, divYieldVector, aMyopicTPlus1Vec);
			
			System.out.println("Fitting Bt+1 at t="+t + "...");
			regBTPlus1Vec = Brandt2005Regression.estimateYValues(excessReturnVector, divYieldVector, bTPlus1Vec);
			regBMyopicTPlus1Vec = Brandt2005Regression.estimateYValues(excessReturnVector, divYieldVector, bMyopicTPlus1Vec);

			// the strategic weights
			optimalWeight[t] = (1+(priceDynamics.riskFree*priceDynamics.dt))*average(regATPlus1Vec) / average(regBTPlus1Vec)/gamma;
System.out.println("At time step "+t+": estimated optimal weight before manipulation: "+optimalWeight[t]);

			optimalWeight[t] = Math.min(1, Math.max(0,  optimalWeight[t] )); //no shortsell and no borrowing constraint
			
			// the myopic weights
			optimalWeight2[t] = (1+(priceDynamics.riskFree*priceDynamics.dt))*average(regAMyopicTPlus1Vec) / average(regBMyopicTPlus1Vec)/gamma; 
			optimalWeight2[t] = Math.min(1, Math.max(0,  optimalWeight2[t] )); //no shortsell and no borrowing constraint
			
			averageExcessReturn[t]=average(regATPlus1Vec); // expected excess return in strategic 
			excessReturnVol[t]=average(regBTPlus1Vec); // expected excess return volatility in strategic 
			averageExcessReturn2[t]=average(regAMyopicTPlus1Vec); // expected excess return in myopic
			excessReturnVol2[t]=average(regBMyopicTPlus1Vec); // expected excess return volatility in myopic
			
			// (TESTING DISPLAY)
System.out.println("At time step "+t+": Estimated Optimal Weight = "+optimalWeight[t] + "(strategic), "+ optimalWeight2[t] +"(myopic)" );

			printComparison(excessReturnVector, optimalWeight[t], priceDynamics.mu0, priceDynamics.sigmaS1, priceDynamics.riskFree, priceDynamics.dt, gamma);
			
			//estimate optimal utility per path per time, i.e. V(m,t), as suggested in p14 Eq. (22)
			if(t==0){
				for(int i=0;i<priceDynamics.pathCount;i++){
					vtVec[i] = utilityFunction(factorVec[i], gamma); 
				}
				
				//estimate V(m,e)
				//regVtVec = LSRegression.estimateYValues(excessReturnVector, divYieldVector, vtVec); 
			}
		}
		
		RealMatrix result = new Array2DRowRealMatrix(priceDynamics.timePeriod-1, 6);
		result.setColumn(0, optimalWeight);
		result.setColumn(1, averageExcessReturn);
		result.setColumn(2, excessReturnVol);
		result.setColumn(3, optimalWeight2);
		result.setColumn(4, averageExcessReturn2);
		result.setColumn(5, excessReturnVol2);
		
		SimpleDateFormat sdf=new SimpleDateFormat("yyyyMMddHHmmss");
		String dateStr=sdf.format(new Date());
		//String outputFilename = "C:\\dev\\MAFS\\mafs6010D\\Project\\strategic-asset-allocation\\optimalWeights"+dateStr+".csv"; 
		//save(outputFilename, result);
		//System.out.println("Result saved to "+outputFilename);
	}
	
	
    public static void save(String fileName, RealMatrix matrix) {
	    try {
	
	        BufferedWriter writer = new BufferedWriter(new FileWriter(fileName));
	        for ( int i = 0; i < matrix.getRowDimension(); i++){  
	        	for(int j=0; j<matrix.getColumnDimension();j++){
	        		writer.write(matrix.getEntry(i, j) + ",");
	        	}
	        	writer.write("\n");
	        }
	        writer.close();
	    } catch(IOException ex) {
	        ex.printStackTrace();
	    }
	}
    
	public static void main(String[] args) throws Exception {
		
		// CRRA utility parameter
		double gamma = 5.0;

		OneAssetSimulationTrial2 t=new OneAssetSimulationTrial2();
		
		//as in Table 1 Panel A
		//VW with Momentum + MR
/*		double phi = 0.15; // degree of momentumness of the stock price process, range from 0 to 1
		double mu0=0.92/100;
		double mu1=0.016;
		double sigmaS1=5.24/100;
		double alpha = 0.011;
		double sigmaX1=-5.77/100;
		double sigmaX2=1.35/100;*/
		
		//EW with Momentum + MR
		/*
		double phi = 0.39; // degree of momentumness of the stock price process, range from 0 to 1
		double mu0=1.15/100;
		double mu1=0.012;
		double sigmaS1=6.28/100;
		double alpha = 0.0094;
		double sigmaX1=-8.43/100;
		double sigmaX2=1.52/100;
		*/
		


		
		double phi = 0.39;
		double mu0=1.15/100;
		double mu1=0.012;
		double sigmaS1=6.28/100;
		double sigmaS2=0.0;
		double alpha = 0.0094;
		double sigmaX1=-8.43/100;
		double sigmaX2=1.52/100;			
		
		
		
		PriceDynamics pd = new PriceDynamics(mu0, mu1, sigmaS1, alpha, sigmaX1, sigmaX2, phi);
		//pd.riskFree=0.03;
		//pd.pathCount=5;
		//pd.riskFreeRate=0.03;
		
		
		pd.simulatePath();
		t.solve(gamma, pd);
		

		System.out.println("Simulation Finished");
	}
	

	


	
	
}
